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ABSTRACT 



The analysis of nonlinear reactor kinetics with Doppler feedback 
is developed by the use of the Poincare-Lindstedt perturbation technique. 
The two types of temperature dependent reactivity coefficients to be 
investigated are 



^ ^ i 

dT T 



and 



~ oc - 
dT ^ 3/2 . 



Each of these reactivity models are analyzed for the reactor with positive 
and negative coefficients of reactivity and infinite slab and finite cylin- 
drical geometries. The Poincare-Lindstedt perturbation solution is con- 
structed in the form of an expansion about the unperturbed state of the 
reactor. Each of the perturbed equilibrium solutions is subsequently 
analyzed for stability. 
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I. THE NONLINEAR TEMPERATURE DEPENDENT REACTOR 



A. GOVERNING EQUATIONS 

In a homogeneous, one-velocity, bare reactor the neutron popu- 
lation is described by the equation (1) 



The first term Eq , (1) represents the time rate of change of the neutron 
flux within a differential control of volume. This term is balanced by the 
net effect of the diffusion of neutrons, the absorption of neutrons, and 
the production of both prompt and delayed neutrons within the volume. 
The concentration of the precursors is described by [1] 



The solution of these equations with temperature dependent parameters 
represents the goal of this work. 

The temporal and spacial coordinates may be nondimensionalized 
by introducing the definitions 



± ’bmA) 

AT ^ t 



D vVrx.t; - 



( 1 ) 






X 
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Equations (1) and (2) become the following 



M 



c 

c. ^ k . 



tit 



AT I, 



( 3 ) 



(4) 



where the precursor concentration and the flux must have a zero value 
on the boundaries of the reactor. The solution is desired for the infinite 
slab reactor and the finite cylindrical reactor. The first of these cases 
is the slab which may be envisioned as a reactor which is infinite in 
two dimensions and bounded by the coordinates +1 and -1 in the third 
dimension. This simple geometry will easily lend itself to analysis, 
but permits the observance of mathematical rigor and the physical inter- 
pretation of the results. The boundary conditions require that the flux 
is zero at the coordinates + 1. The second case is that described by 
a finite cylinder which has boundaries on the longitudinal axis at + 1 

and at a nondimensional radial distance of r . Therefore, the flux is 

o 

required to be zero at the extremes of the longitudinal coordinate, 

Z = + 1 , and at the radius r = r . 

~ o 

B. DOPPLER EFFECT 



The lumped parameter k is generally referred to as the neutron 
multiplication factor. In elementary reactor analysis the multiplication 
factor is evaluated as a constant given by the product [1] 
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^ Pt 



However, in reactors that operate at high temperatures the multiplication 
factor becomes temperature dependent as will subsequently be shown. 

One of the major considerations to be made when accounting for 
the temperature dependence of the multiplication factor is that described 
by the nuclear Doppler effect. Physically, this phenomenon results from 
the temperature broadening of the capture cross-sections of the fuel 
nuclei!. Thus the Doppler effect may be envisioned by a temperature 
dependent competition for neutrons by capture and fission processes. 

The quantitative description of the Doppler effect is achieved by making 
the following assumptions: 

a. The fuel density as a function of velocity is given by a 
Maxwellian distribution function. 

b. The capture cross-sections near an isolated resonance is 
given by the Breit-Wigner formula . 

It may be shown that the average radiative cross-section, CT ^ ^ 

and scattering cross-section, CT $ , are described as a function of 

energy and temperature given by [1] 




and 
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where 






and 



C ^ > X } are tabulated functions of 



temperature , 



A 

a 



2 

r 



(E -e.) 







The effect of C3* y and O* $ is most conveniently incorporated 
into one quantity by introducing the effective resonance integral, I, 
defined by [2] 



r 




EFFCCTivt 




f ^ dE 

J Olga's E 



where cr» is the additional scattering cross-section per fuel atom. 
It is usually more convenient to analyze the temperature dependence 
through the J function which is related to the resonance integral as 
follows [1] 







J(X'P) 


where 




V 




■ a'a cr, 


T* 

' Y 



and the subscript M refers to the moderator. 

If attention is returned to the concept of the multiplication factor, 
its dependence upon temperature may be investigated [1] 

Since 

i -p 
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and 



_>p Q/)Cp (I) > 

the Doppler effect on k is studied by considering the temperature effect 
on . However, the preceding analysis showed that 

Om r* V 



I = 



j (tf) 



and it may be simply concluded that 






oc 






J^b|3)] 



If the derivative with respect to temperature is taken, it is obvious that 



1 ci ^ d J 

o' — oc • 

R dr (±T 



(5) 



However, under most operating conditions the value of ^ is near 
unity. Thus for small changes in ^ the following relation is valid. 



d i ^ d J 

dT dT 



(5a) 



The left hand side of Eg. (5) is referred to as the temperature coefficient 

of reactivity and provides a functional dependence upon the 3" function 

and, in turn, upon the changes in temperature of the reactor. 

In a paper by A. Reichel the temperature dependence of the 

function was investigated [2]. The results were based upon evidence 

238 

gathered from experimentation with U and from analytic investigations 
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into the nature of the 



J 



function. Reichel concluded that physically 



realizable states of the reactor must be located in the range 
.02 ^ ^ ^ ’(o and . OOOi < /3 < 10 • 



The results are most simply presented by Figure 1. The variable 
described by K where (l 

( 2 ) 

The shaded region represents states which are physically realizable for 



/3 



IS 



238 

U . From the analytic investigations it was observed that values of 
^ at the extremes of its range coincided with conditions which 
required that the 3 " function be a constant.' Therefore, the tempera- 
ture dependence of is destroyed and, in turn, the Doppler effect 

and temperature dependence of the coefficient of reactivity is nonexistent. 
In other regions of Figure 1 the dependence of X was found to vary 

with temperature raised to a constant power. These regions are indicated 

cons t 

on Figure 1 by the lines labeled with T . Thus for any particular state 
of operation defined by values of and ^ the dependence of 

upon temperature may be defined. If Eq. (6) is evoked, this T de- 
pendence is transferred to the governing equations, Eqs. (3) and (4), 
via the temperature coefficient of reactivity. 

As may be noted from Figure 1 the variance of the X function 
with temperature may take on many different forms. However, in an 
analysis presented by Hummel and Okrent [3] the X function was 
found to vary betv/een the functions 

~'/z 



J ^ 



J" OC T 
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jr (tip) 



FUNCTION TABLE 



I 




Figure 1 
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I 



These results are used to determine the functional dependence of 

with temperature. The results are that 

-1 






and 



dt 

dT 

dT 



oc 



oc 



T 



T 






(6) 



1 



respectively. In the remainder of this paper the effect of these types 
of temperature feedback are analyzed. 

The proportionalities in Eq. (6) are removed by the addition of a 

d 

constant. The values of determined by experiments per- 

formed on the reactor in question. Some representative values are given 
for various reactors . 



TABLE 1 [4] 



Reactor 



Nuetron 

Lifetime 



Temperature 



Coefficient 



Fermi 

EBR II 

Hallam 

Peach Bottom 

EGCR 

EBOR 



14 1 ^- sec 

68 yU sec 

3 7 m sec 
23 m sec 

4 0 kv\ sec 

22 ^ sec 



-3.4 X 10 
-4.2 X 10 
+ .3 X 10 
-2.6 X 10 
-4.5 X 10 
-1 . 5 X 10 



-5 

-5 

-5 

-5 

-5 

-5 



h k/k per °C 
^ k/k per °C 
^ k/k per °F 
A k/k per °C 
A k/k per °C 
A k/k per °C 



^In the subsequent analysis the need to differentiate between the 
feedback models will arise. The distinction will be made by referring 
to the "T"l model" or "T“3/2 model." 
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model reactors the temperature dependence of 



-1 -3/2 

For the T and the T 



a 

T 



the multiplication factor may be described by the equalities 

_ 

d T ~ 

and ^ (7) 

d R _ _ ^ 

dT " J’A 

respectively. It may be noted that the negative sign is incorporated 
into above equation. This convention is generally observed because of 
the greater frequency of reactors with negative temperature coefficients . 
In Eq. (7) the units of CX, and 0. are defined such that ^VdT 
has the units of ^ k/k per °C. 



C. ASSUMPTIONS 



Before the analysis of the reactor may be undertaken two further 
assumptions must be made concerning the dependence of the neutron 
flux with temperature and the delayed neutron concentration with time. 
The first assumption states that the reactor temperature rises instantan- 
eously with the flux (this is the usual prompt feedback model) . Equation 
(7) becomes 




a 



a‘f> 

and 

d ^ 



g. 



( 8 ) 



The second assumption states that the delayed neutron concentration 
has delayed neutron concentration has reached its steady-state value. 
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Delayed Neutron 
Group 



TABLE 2 [1] 



Half life (sec) 



Br 



87 



54.5 



I 



137 



24.4 



Br 



88 



16.3 



I 



138 



6.3 



Br 



89 



4.4 




6 



The justification of this assumption lies in the comparison of the neutron 
life time in Table 1 with the precursor half-life in Table 2. It is evident 
that a number of neutron generations will have progressed before the pre- 
cursor concentration has been changed. 

The assumption of steady-state precursor concentration reduces 
the neutronic equation to the case without delayed neutrons. This work, 
therefore, deals with the case where either no delayed neutrons exist 
in the reactor or the delayed neutron concentration instantaneously follows 
the flux distribution. The latter case is realistic when considering the 
reactor response to small perturbations. 
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II. ANALYSIS OF THE REACTOR WITH A 
TEMPERATURE COEFFICIENT 



^ - J- (T ^ Model) 

dT T 

The solution of the nonlinear equations, Eqs. (3), (4) and (8), 
governing the reactor kinetics is accomplished by means of the Poincare- 
Lindstedt perturbation technique [5] . This method offers a solution to 
the neutron flux distribution for small perturbations in the flux intro- 
duced by the temperature feedback effect. The perturbed flux distribution, 
referred to as is described by a Taylor series expansion 

about the initial state where € is the amount of perturbation 

from ^ . The perturbation technique yields multiple equilib- 
rium solutions for and associated eigenvalues for each 

of these states. The presence of more than one equilibrium state, known 
as the bifurcation phenomenon, requires the investigation of the stability 
of each of the possible states of flux. Stability considerations are made 
by imposing an initial deviation on an equilibrium state and examining 
the resulting temporal behavior of this state. 

A. EQUILIBRIUM STATES FOR A SLAB REACTOR 

Before the general solution to Eqs. (3), (4) and (8) is undertaken 
the equilibrium states for the slab reactor is investigated. For the case 
where the flux does not vary with time (i.e. ^ ) and the 
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reactor kinetics are not dependent upon delayed neutrons the governing 



equations reduce to the following boundary-value problem in the special 
domain D[-l , 1] : 



^ C'k') = O 



l/n D 
on S O 



d i(V) ^ 

( 8 ) 

k(^o) = go • 



Equation (8) may be solved to give 







a Ln. 



^oCx) 



(9) 



This result may be substituted into the boundary-value problem and the 
result is 






- a 



< ^(x) ' 






I ^(x) = o in D 
(10) 

on S b 



where 





1 



The solution to Eq. (10) is constructed about a flux distribution 
and eigenvalue ^ o which will be determined. 
Physically, Ao is this eigenvalue corresponding to the initial 
material buckling of the reactor before perturbation by feedback. After 
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the system is perturbed the flux and the material buckling establish a 



new state which may be described by a Taylor series expansion about the 
initial state. If the amount of perturbation is measured by € , the 

resultant reactor state is 

= VoM + ^ -f 

'XCCz) = 



Thus for small perturbations of the system the values of 
and X<t) may adequately be described by truncating the series after 
the third term. The unperturbed flux, % ('X) , and eigenvalues, 

Xo > defined as the solution to the linear problem 



+ ( Ao- o.) = O 

^ ( X ) = o 



In P 

( 12 ) 
(5 P 



Since 




and a are constants, the solution is 



and 



Cy.) :=■ Aa^ Cos (X "X 

Ao- a = 



(13) 



where /n = 3, 5 For the purpose of simplifying the subsequent 

analysis the shape function for the slab reactor is defined as 





y) 



(14) 
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The remaining terms of the expansion given by Eq. (11) are determined 
by differentiating Eq. (10) with respect to € and taking the limit as 

^ tends to zero. The values of ^ and A are then determined. 

* * . . 

By repeating this process the values of ^ and ^ are evaluated. 
Thus, differentiating Eq. (10) once yields 



+ A¥’ + - a4> ^ o. aa 



The limit as ^ tends to zero requires that 

C o 



J?/ 



m. 



^ ) — A 



^ ^ o 

As ^ tends to zero, Eq. (15) reduces to 



vVc«) a) Vtx) = in i (16) 



- Q 



on 



Sl> 



If Eq. (16) is to have a nontrivial solution, the solvability condition 
(provided in Appendix A) may be evoked. Thus 




(-A tc-x)) 



dy. 



O 



where ^ C'X') is the homogeneous solution to Eq. (16). The 

f-l c 'n ^ 

solution for Cx) is easily seen to be ^ ("k) — ///j . 
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When the functions 



are substituted 






and 






into the solvability condition the result is 



-A \ ch( - o . 

^ t> 

Since the integral is nonzero, it follows that x=o . This result 

m 

is substituted into Eq. (16) and the solution for (■PCX) is clearly 
seen to be ^ (x) = e^ . 



Thus the linear terms of the expansions 



are 



— S/n ( “X) and 

A = o . 



The parameter ^ is chosen to represent the average magnitude of 
the linear term in the expansion of Eq. (11), weighted by the 
distribution, and measured in units of the maximum initial flux. Thus, 

e dfx . (17) 

Differentiating Eq. (17) with respect to ^ yields 

f\^ - \ ^C'x,k) c/x 

Jn 



where 



- cpcx) + ^ Cpix) + 
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By letting £ 



approach zero 





— 




CM 


d')( 






Jt> 




thus A/yj - 


Bai 


. The second terms of the expansions are 


1) 


/^At 


'tr 


and 


11 

0 

• 


The third term of 


<^P(X,(r) and 


^ C4 ) are similarly 



determined by differentiating Eq. (16) with respect to £ . The result 

is 

+ Mo + 2M<^'> ^cxjO 

- a - a ^cx,i) = o i" 

^ C\e) - o on5l>. 

Since A=o , the result of taking the limit as tends to zero 

is 

S7^^(X) + (M~ 0.) ‘4>iX) = (19) 

* ^ ^ - 

' 4 > Cx) = o 



- in J> 



on . 



The solvability condition requires that 




d ^ (X) 



1 






^ (x) 



dx 



O 
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• • H 
V (X) 



where ^ (X) is the homogeneous solution to Eq. (19). This 

* ‘ H H / 

solution is obviously seen to be = C X, 

When the various functions are substituted into the solvability condition 
the following becomes evident 



X 



c/x = 

• I 

and X = a . This result is replaced in Eq. (19) and it is found 
that the right hand side identically becomes zero. Therefore, the solution 

for the third terms of the expansions is ^ (x) - C/n and (2 0) 

1 1 

A = Ol . The equilibrium solution for ^P(Xj^) is obtained 
by placing the results of Eqs . (13), (18), and (20) into the series of 
Eq. (11). Thus, 



= /l^ (i 2 ^ (X) (21) 

and 

\U) - <X = ('^^2 )^ + 2 CL + • 

where = 1 , 3 , 5 . . . . The value of from Eq . (2 1) may 

be evaluated if the constants and are defined. Since the reactor 
is assumed to operate at its fundamental mode before temperature feed- 
back becomes effective (later it will be shown that the other modes are 
unstable) , the value of may be determined by the power requirements 
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of the reactor. Thus, 



P 




d-x 



where P is the linear reactor power level and is a conversion 

factor. The value is undetermined and will remain an unknown 

z 

of order € throughout the analysis . For each mode, /VI , 
as defined in Eq. (21), describes two possible states in the reactor 
as depicted in Figure 2 for the fundamental mode. The state of larger 
flux is determined by a positive value of € while the smaller state 
is described by a negative value of £ . Associated with each state 

of the reactor is an eigenvalue, “* CX. , given by Eq. (21) and 

plotted against G in Figure 3. The values of ACO - ^ for the 
case of negative reactivity are drawn as a solid parabola which opens 
to the right. It may then be noted that the equilibrium state may only 
exist for values of A(0- a which are greater than the value of 
Xo ' initial material buckling of the unperturbed system. For 
the reactor with positive reactivity the converse is apparent. The cor- 
responding curve for positive reactivity in Figure 3 is that drawn with 
a dashed line. These eigenvalues are found to be always less than the 
value of Xo - It n^^y also be noted that the eigenvalues of the 
system are A(e) - CL where the term ( - a ) describes the 

amount the material buckling changes so as to accommodate the temper- 
ature feedback . 
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FLUX DISTRIBUTION FOR SLAB REACTOR 




Figure 2 
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REACTOR BUCKLING 




- - - - negative reactivity 
positive reactivity 
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B. STABILITY ANALYSIS OF THE SLAB REACTOR 

The stability criterion for a slab reactor is developed by considering 
the temporal effect of an initial deviation on . Under 

these conditions the governing equations, Eqs . (3), (4) and (9) reduce 
to 

= yV(X;Tr) + ^ 

"h 'X. 

^C‘f) = - a JK 

V(^o) = = o orx£^. 



These equations may be compacted into the more convenient form 

- a 

■= o 



c)^CX,T) ^2.^^ 

= V ^CX,X) -h 



Nl't 



9 (va) 



oC^), 

In ^ (22) 

on ^ ^ 



where 



A = - 1 . 



The time dependent nature of the neutron flux is investigated 
by prescribing an arbitrary initial deviation, , on the equilibrium 

flux and then observing the growth of the.deviation . An unstable state 
of flux is that which will result in an unbounded growth of the deviation. 
A stable system will behave in a neutrally stable manner which implies 
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that the deviation tends to a new equilibrium state or in an asymptotically 
stable manner which implies that the deviation decays to zero. The time 
dependent flux is described by a Taylor series expansion about the 
equilibrium state q>(x,o , defined by Eq. (21), 

I 1 (23) 

such that 

iim. . 

o 



Equation (22) is differentiated with respect to ^ . Therefore, 









OYk S D 



where 



The limit as 



tends to zero is taken and the result is 



' b C y, e, t ) 
*2> ^ 

- a 













CX (24) 

lV\ D 



r) - o 



Ojrv S J> . 



28 



If a solution for 't ) is chosen to be 

= u.cx,e.) &xp(oCu)'t:) as) 

then the stability of the system is analyzed by inquiring whether the 
value of oUi) is negative or positive. A form for CO 

is given by the series 

CC;) - o(o "*■ £ ^ 2 ^ ^ 

The value of ol ( ^) is determined by the perturbation tech- 
nique previously used. In this vein, Eq. (25) is substituted into Eq. (24) 
and the resulting equation is 

+(>c^)-a - olC(:) - ( 27 ) 

0 \ 

(X )Uy\ Cx) / IXCX^t) •=■ o m x> 

(-A. CXy^) ~ O Oh Sl> , 

Setting £ equal to zero 

/ X (28) 

V U Cx^t) - ( Aoai - CX - oLo) UCx,o) = O m D 

U ( O) •= O orv 6l> , 

Since the values of X o/h ^ ) have been prescribed by Eq. (13), 
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Eq. (2 8) may have a nontrivial solution only if (proof is given in 
Appendix C) 

A. 0/yr\ Q. A OrA ““ Q. OL o • 

Thus the value of oLo is defined as 

O/Vl Av\ ““ ^ A O/vn ^ A O/Vvv “* (2 9) 



- ( y2 )* 


( ' Arv? ) 


for odd intergers of /V^ and /W^ where ywi is the mode of the 

deviation. If the fundamental mode of flux distribution, At i , 

considered, it is evident that Eg. (2 9) reduces to 


dcxf^ - ( "% ) 


( 1 “■ ) 


Thus the values of o(,o\ /m ^re 




d. oi ~ O 


/m — 1 


- - 2tj 


/m = 3 , 5, 7, • • 
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For values of ^ near zero the first term of the series for otco 
indicates that the fundamental mode is not unstable for any value /m . 
However, the neutral stability predicted by = i obliges an investi- 
gation into the higher order terms of Eg. (26). Some conclusive results 
may be obtained by considering the linear term of ( O . 

For higher modes of the reactor (i.e., /r\^i ) Eg. (29) indicates 

that 

~ O /vw = A\ 

— AVV AC 

Since a positive value of o<(e) may occur, it is concluded that the 
higher modes of the reactor are unstable. For a physically realistic 
situation, however, only the fundamental mode exists in the reactor 
before the perturbation. 

The solution for U O) is obtained by substituting the 

value of Eg. (29) into Eg. (2 8). The result is that 



S7 UCOC^O) + ( Aoatx " L^CX,0) = O in D (30) 

UC*K,0) - O on 6£>. 



It is then obvious that 



La. c "X ) o ) \X 
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where the index (TA designates the modes of the deviation. The value 
frr\ = 1 refers to the fundamental mode of the deviation while > 1 
refers to the higher modes. 

The linear term of the expansion for oca) , defined in 
Eg. (2 6), is derived by differentiating Eq. (2 7) with respect to € ; 

thus 



V" ucx,6:) + (ac€.'> - a - oL(^) 



(31) 



- a 

Ca) ) 



- a 



UlO<,^) = O 

The limit as ^ tends to zero is taken and 



I n D 
O Kk S b • 



u.(x,o) + ( X 



onw 



- a 



^OoTt/vn ^ 



(32) 

U*('X,0) = 



i-i 



Co) + cJdco) + a 






9iCy() 



U(X^O) in D 

UlCX^o) = O ov\ 

By noting that A “ O the solvability condition for Eq. (32) 

requires that 




od + CL 



CPCA) 



\ . H 

J Uw (x,0) UL iXjO) 



dlx 



= O 



. H 

where U ( X^o) is the homogeneous solution to Eq. (32). When 
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the value of / from Eq. (29), is employed the function 

• W 

U Cx,o) is determined to be 

U Cx,0) = U 'V'/yv, , 

The value of <X is defined by substituting the functions 
^(x) (from Eqs . (13) and (18)), LXCX^O) , and (JL (XjO) 
into the solvability condition and the result is 

^ ~ - Ol . (33) 

This result and Eq. (29) determine the linear approximation to 
^ ) » Eq. (26). Thus, 

oLi^) = (■%) (pr{ - rm) ~ € <X 

It now becomes evident that the question posed concerning the stability 
of the fundamental mode (/yv-i ) for sufficiently small C is 
answered by the linear approximation for C>{ ( ^) . Since values of 

/yy\ > i give very large, negative values of ^ (€) and thereby 
assure the stability of these modes, the concern is centered about 
/7V| ■= /yi ” i • For this case 

0 ^ ( ^) = ""Ccl 

and the stability of the fundamental mode will depend upon the sign 
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of the linear term of 0^(0 . The various cases to be analyzed are 

the following: 

a) negative reactivity 

oUi) < O for positive € 

> Q for negative ^ 

These equilibrium states of flux are depicted in Figure 2 and the corres- 
ponding eigenvalues are shown by the solid curve in Figure 3. A positive 
value of C describes an equilibrium state of flux which is larger 
than the initial flux. This state, Eq. (21), exists for an eigenvalue, 

- a , greater than the initial material buckling of the unper- 
turbed reactor, Xo ' observed to be stable. A negative 

value of £ indicates an unstable equilibrium state in the reactor. 
However, Eq. (21) indicates that the eigenvalues must be greater than 
the value of Xo for any value of £. . Therefore, a situation in 

the reactor which permits a decrease in flux (i.e. , negative £. ) with 

an associated increase of the reactor buckling is physically unacceptable. 

b) positive reactivity (evaluated by changing the sign of all terms con- 
taining a ) 

ck(C) > O for positive ^ 

< o for negative G 
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These states of flux are depicted in Figure 2 and eigenvalues shown by 
the dashed line in Figure 3 are analyzed in a similar manner. The stable 
state of flux is that described by a negative value of £ . From Eq. (21), 

noting that the sign of CL must be changed to accommodate positive re- 
activity, the fundamental eigenvalue is less than that of the initial 
material buckling, X© . Therefore, a decrease in flux (i.e., negative 
£ ) is compatible with 'a decrease in the eigenvalues. However, for 

the case of positive £ the equilibrium state is unstable. Again it 
is noted that the increase of flux associated with a decrease in the 
material buckling is uncompatible and this situation is physically un- 
realistic . 

The solution to Eq. (32) can be completed by first substituting 
Eq. (33) into Eq . (32). The right hand side of this equation becomes: 



RHS + a UC'X,o) . 

However, from Eqs . (13), (18), and (30); 







respectively. Thus, 



RHS 




O . 



Equation (32) correspondingly reduces to the form 
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LU^O) + ( Aom- o. - = O I’ 



UCX, O) = o 

0 

Thus, the solution to LA-CX^O) is 



On 4 D . 



. H 



Cx,o> = . 



(34) 



U(X,0) = LL _,vv, ,AK 

For larger arguments of ^ the thircd term of the expansion of 
0^(0 , Eg. (2 6), may become important. This term may be evalu- 

ate(i by (differentiating Eq. (31) with respect to ^ an(d then taking the 
limit as ^ ten(ds to zero. The differentiation yields 

v" + ( Xt) - a - o2(6) 

" ^ ) u. + 2 ( 

- oiU) - 



a ) ui.cx,ii) 



+ - ck U) - a 



-\r (X. 



(4> C'K.fz) 

<f"cx 



^ c A; 6 j 

) U C 7 <j^) = O ^ 



U.(X^6) = 0 ovn^T*. 

When the limit as ^ tends to zero is taken and the various known 
functions are substituted into the resulting equation, the following 
becomes evident. 
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^^7 LXCX^O^ ^^o/v\ ^ ^Om/vr^ O-CXjO) 

[u - a (2 - fc)] 



(3 5) 



uca,o) = o 

The solvability condition requires that 



OK 



\n D 



S D . 



I [^ - a (2 - aZ)] d')( =0 



H 



where U CX,0) is the homogeneous solution to the above equation. 

'• H H .. H » 

The nontrivial solution for LX C%0) is UL Cx,0) — »S»»v • 



Thus , 



[c^- a ( 2 - §T^ 1 / d'X = 



/\ /M '^/>n * 



o 



and since the integral 



[ ild'x-i. 

^ = a'( 2 -fe) 



(36) 



The expansion of d. ie) , Eq. (26), may be written in the more 
complete form by the substitutions of Eqs . (29), (33), and (3 6). Thus 

0( (e) (rn-orr!') - £ 0. 



+ 



C, 



(2- s: 



)+■ 



( 37 ) 



For larger arguments of ^ the stability of the system is contingent 
upon the negativity of oi (e) , defined by Eq . (37). For the 
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I-I 



fundamental mode ( /Y\= i ) Eq. (37) reduces to 

oUi) =(yif (l~ryn) - £(X (2- fe,) + - 

The value of /m of particular interest is that given by fW\ ~ 1 since 
values of /Wl > 1 provide a large, negative first term in Eq. (38). 

Thus, the stability criterion is 

+ %,) < O 

or 

- £ a f 1 - £ + i € ^>^, ) < O (39 

By inspection of Eq. (39) the quantity within the parenthesis may be 
seen to always be positive for small arguments of ^ . Thus, the 

stability is governed by the coefficient ( •" £ Q ) as has previously 
been considered. However, for cases of larger arguments of ^ the 
inequality shown above may be changed. Thus, to the order of t 
the condition for stability is 

1 

^ 1 - i Ya, ■ 

If this were violated, such large values of ^ would render the re- 

« • 

actor unstable. Thus, the undetermined peak amplitude of 9 (ax') , 

C^, as compared to the peak amplitude of the unperturbed flux, A^, 
is important when considering the stability regime. 
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C. EQUILIBRIUM STATES FOR A FINITE CYLINDRICAL REACTOR 



The finite cylinder represents a realistic geometric configuration 
for a reactor. The physical and mathematical considerations which have 
been deduced for the slab reactor are directly applicable to the cylindrical 
reactor. In this section, the equilibrium states of flux, associated 
eigenvalues, and the stability of these states which have undergone 
feedback perturbation from some initial state will be studied for cylind- 
rical geometry with azimuthal symmetry. The flux distribution must 
depend upon two special coordinates identified by P , the radial 
coordinate, and ? , the longitudinal coordinate. 

The equilibrium solution for the reactor is developed from the govern- 



ing equations under the conditions that 

The governing equations reduce to the form given by Eq . (10) which is 
rewritten below, 



The feedback-perturbed solution to Eq. (41) is evaluated by a 
Taylor series expansion about the initial state. If the amount of per- 
turbation is ^ , which will subsequently be defined, the resultant 

equilibrium state is 






o 



on , 
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( 42 ) 



-i- + z^ ?\ "^ ' ' ‘ ^ 

The initial state , and Ao /is defined as the solution 

to 



V^f„(f.2) + (A.- a) = O m s M3) 

( f"; ?: ) " O OK 5 - 

Equation (43) is solved by assuming a separable solution for 
given by 



— R (r) ^ C^) , 



(44) 



Substituting Eq . (44) into Eq . (43) yields the following: 

J- d / d Rcr )\ 
Rw d ■*■ ?'^2) Tc*r('^ av- J 

+ ( Ao- a) Rcr) = O 



which may be separated into the equations 



/ c^R(r)\ 2 o 

dTr^ J + V* Ui R(r) 



d r 



O 






R(r) — O 

^ 2 C2^) = O 

Z^tL) O 



Zci ) 

d e’- ^ 



c5v-» ^ 

I v\ 

0\rv S I> - 
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— Z Z 

Co + Co — Ao *" O. . 



The solutions to these equations are the following: 





Rcri 


zz 


m do %■<, ) > 








ZZ 


Cos (^J /2 i) , 




and 


Ac - a 


= + (^7/2 )' 




where 


fJ-/^ are 


the I 


■oots of the zeroth order Bessel function. The 




solution to Eq. (43) is 


; » 










RCr) ’S C^) 










A/n Jo (j^r^V^c) CoS ^ 


) 


and 


Ao- 


Ql 




(45) 


where 




=s 


i. A • • • . 




For the 


purpose oi 


[ simplifying the subsequent analysis the shape function 




is defined for the cylinder as 






i'm 




J. Cos h). 


(4 6) 


Thus , 
















= Aa, . 


(4 7) 
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The remaining terms of the expansion of and 

Xa) , given by Eq. (42), are determined by differentiating Eq. (41) 
with respect to ^ and then taking the limit as ^ tends to zero. 
These steps lead to the equation 



-f (;\o-a) == - A m 

^ C^,^) - O OK (48) 



If is to have a nontrivial solution, the solvability con- 

dition (provided in Appendix A) may be envoked . Thus, 



f • * H 

I - ^ 2jj‘Y' dr d *2: = 



• y/ 

where (r,i) 



is the homogeneous solution to Eq . (48). 
solution has a similar shape function as Eq. (43). Thus, 



o 



This 



(49) 









(50) 



/O r* 

The fundtions ^Cr, and ^ Cf~j a , Eqs. (47), and (50), 
respectively, are substituted into Eq. (49); hence, 

/^ - X A« 'I B,:! 2j r dr d H = O 

or 

- % Am B/yv 2j'j f dr ’== O . 
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\ 



I 



i 




Since the integral 



r dr dz = j 1^0 i^„) 

^ O j 

the solvability condition may be satisfied only if A “ O . This 
result may be substituted into Eq. (48) and the value is 

the solution to 

K7"<p(r;i) + (Ao-a)V^r,j; =0 J> 

» 

= o . 

Hence, 

4^(r,i) = Cf>(r,i) = , (50a) 



The parameter ^ is chosen to represent the average magnitude 



of the linear term in the expansion of Eq. (42), weighted by the function 



3T Y'J’ 

Thus , 

1 



and measured in units of maximum initial flux. 



f = A 



/V\ 



JJ y^o Ji ^jZrn) -* 



(51) 

Zjv'dr d"£ 



Differentiating Eq. (51) with respect to ^ and letting ^ approach 
zero yields 
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f\ 







2 t t r dr dn 




Therefore the linear terms of the expansions of and 

^ C ^ ) are 



(53) 



and 



A = o. 



The third terms of the expansion of and A 

are determined by differentiating Eq. (41) with respect to ^ twice 
and letting ^ approach zero 

+ (A»*a) - 

- o 

The solvability condition requires that 



(54) 

m I> 

CsV-N. ^ E> . 



/ / (p (Ci) V /A \ 

J (Cl ^ ^ (rji) 2yr dr c\^ ^ O 



/; H 



where ^ (!<■/ ^) is the homogeneous solution to Eq. (54). The 

' * U 

solution for is similar to that obtained from Eq. (43), 



thus 
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, determined from Eqs . 



cf r.o = cl . 

The functions and 

** 

(47) and (53) respectively, and Cp are substituted into the 

solvability condition. Thus 

l(aA^t- drds=0 

{a-\)}^flrdrdi-0. 

Therefore, 

• • 

}\ - (X . (55a) 



If this result is placed into Eq. (54), the right hand side will identically 

i» * 

be equal to zero. Thus the solution for is similar to that 

obtained for Eq. (43) and 








(55b) 



The various results for the terms of the expansion of Eq . (42) are 



presented as follows: 



(l + C + i c' 

7 7 (56) 
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It is apparent that for any mode of flux there will be two possible 
equilibrium states defined by the negative and positive values of ^ 

For each of these states a unique eigenvalue is present. For negative 
reactivity ( ) must be greater than Xo positive re- 

activity the opposite must be apparent. It may be noted that the eigen- 
values are defined by - a . Where the term - a denotes 

the amount of change of the material buckling to accommodate the tem- 
perature feedback. 



D. STABILITY ANALYSIS OF THE FINITE CYLINDRICAL REACTOR 



The stability criterion for a cylindrical reactor is developed by 
considering the temporal effect of a deviation on the equilibrium states. 
The governing equation has the form 



c) 

c) 't 




cr,i) . 







(57) 






The temporal nature of the flux is investigated by imposing an 
initial deviation, 9^ , on the equilibrium states defined by Eq. (56) 

and then observing the growth of the deviation. An unbounded increase 
will indicate instability; neutral stability implies that the deviation 
tends to a new equilibrium state; and asymptotic stability is assured 
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by the decay of the deviation to zero. The deviation is described by a 
Taylor series expansion about the equilibrium state , de- 

fined by Eq. (56): 

+ -VJ ^ Cr,i,i,'i) (58) 

such that 

iiKvi _ 

'Vj —>o 



Equation (57) is differentiated with respect to and the limit as 

tends to zero is taken. The result of this process is 

+(>(0 -a-aS^ 



(58a) 

t) iw 



% =0 it) . 

If a solution for is chosen to be 

Lii'Ciji) QMp ( o{(^) 



then the stability of the system is determined by the sign of 0^(6) . 



A form for 



oC(0 



is given by the series expansion 



0{(^) — O^o 2. & ^ ■*“ 



( 59 ) 



The value of (^(0 is determined by the perturbation techniqu( 
previously used. The function given by Eq. (58b) is substituted into 
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Eq. (58a) and the resulting relation is 



V' u + ( > (4^ - a - CO - 0. J!/y^ 






(60) 

U =0 






= o 



OV\ 






The value of (X® may be determined by taking the limit of Eq. (60) 
as ^ tends to zero. Thus, 



(61) 



y^UCir; (X ~ oi^ ) U cr, 1) = O in X> 



U Cir; z ) = O 



On •i'ti 



A nontrivial solution to this equation is obtained only if (see Appendix 
C for proof of this) 

- 0. = X.r. - a- c^o . 



Thus oCa is given by 

^^o/nrvn ~ {Ao/w"” ^ A o/Tvi ■“ 

Some remarks may be directed toward the stability of the system via 
Eq. (62). As was apparent for the slab reactor the fundamental mode 
( /VI =1 ) may be evaluated as follows: 

O(o,^ = O /m = 1 (62a) 

< ■ 2 TJ- ^ A-n > 1 . 
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Thus the higher modes of deviation (/W^l ) which are imposed upon 
the equilibrium flux will vanish with time. However, the value of ~ 1 , 

provides O(oirn\ ^ and the system may be thought of as neutrally 

* 

stable. Thus it is necessary to inquire about the effect of oC upon 
the system for the value of 1 . 

The higher modes of the reactor ( /Vi > J. ) are analyzed as follows: 



ck (fe) < 




/W) 


> /y\ 






o 


/w\ 


- m. 


(62b) 


> 




kVVV 


< /V\- . 





Since the condition provided by /Wl ^ /A may be realized for modes 
of /yi ^ 1, , values of o((^)> o are present and the higher modes 

are unstable. However, the reactor is assumed to be operating at the 
fundamental mode before the feedback is introduced. Equation (62) for 
Oiomrm. iT^ay be used in Eq . (61) to provide a result similar to Eq . (43). 
Thus , 

'm- («) 

(See Appendix C for the proof of this.) 

The linear term of the expansion of (Xfe) (Eq . (59)) is evaluated 
by differentiating Eq. (60) with respect to ^ and taking the limit as 
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^ approaches zero. The result is 



+ (> 0 ^ - a - (XorrvrvK) (64) 

= Od + a ucr,i) irt D 

uCr^ 2 ;) = o GA c5:i) . 

By noting that has previously been determined to be zero the 

solvability condition requires that 



J r / • \ • ^ 

' {cK + a ucr,*) = o 

7\ 



(64 b) 



. u 

where U is the homogeneous solution to Eq. (64). By the 

same reasoning used for the solution of LL the solution for 

. H 

is 

■ ^ - TT AP 

LL Cf, s ) LA. i/m. 



When this result and those provided by Eqs. (47), (53), and (63) are 

« 

used to evaluate Eq. (64b) a value for is derived . Thus , 



If. 



oC + a 



4^ 





dz = 0 




O . 
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Therefore, 



(^ = "■ 0l . (65) 



This result may be substituted into the expansion for o( (^) and 
further analysis of the stability is made. Thus, 



o( (^) = o( 



om rm 






( 66 ) 



The question pertaining to the stability of the fundamental mode (m ~ l') 
is resolved for small ^ . For the case where /yy\ = 1 the value of 

0( oil is zero and 

o( (i) = - 6 cx . 



The analysis of the various cases of negative and positive feedback is 
similar to that developed for the slab reactor. The conditions of stability 
for the fundamental mode are briefly stated as follows: 

a) The stable equilibrium state for a reactor with negative reactivity 

has a positive value of ^ . The other equilibrium state is unrealistic. 

b) The stable equilibrium state for a reactor with positive reactivity has 

a negative value of ^ . The other equilibrium state is unrealistic. 

If the value of cC is substituted into Eq. (64), the right hand side 
will become zero. Therefore the solution for u is the 

homogeneous solution to Eq. (64) and 

yjL. ( r, ) =: XX rvy\ ^ /nv - (66a) 
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For larger values of ^ the quadratic term of may be 

necessary. To obtain the desired results Eq. (60) is differentiated twice 
with respect to ^ and the limit as ^ tends to zero is evaluated. 

If all of the known functions are employed the result is 

+ (XoM ~ <3. - =(67) 

(ck - a (2 -^Va^)) in 

LA-C'f', — O . 



The solvability condition requires that 



j[^ - a (2 

# # H 

where IX C 




“ ^ )] U ir,^) Ijr dr di^ 

is the homogeneous solution to Eq. (67). Therefore 
. It is obvious that the solvability criterion is met 



only if 



or 



oC-a(2"/0 =0 

ck = a ^2 - Va.') . 



( 68 ) 



This result is substituted into the expansion for ) , Eq. (59), 

and the following becomes apparent, 
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, X Z , V 9 

0 (( e ) (h-1 -m^) 

-€a +j^"Ql(2- Va„) +'••. 



Since the fundamental mode of operation is generally of greatest interest, 
Eq. (69) is rewritten for /^l “ 1 . 



1 ^ 

(X — ( To ) (/^i ~~ /^/m } ( 1 ~ /yy\ ) 

-€a[l-l€(2 - %,)] + •■•_ 

For /VV[ > i the first two terms dominate and cK Ce) < 0 . For 
/m = 1 the stability may be determined by evaluating 

o((e) = -€aj” 1- z^- {2 - ^7^^^ . 



These results are similar to those obtained for the slab reactor and 
stability is granted under the conditions established for Eq. (37). 
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J. 







III. ANALYSIS OF THE REACTOR WITH A 
TEMPERATURE COEFFICIENT 



^ = - a T MODEL) 

dT 

The solution of the nonlinear equations (3) , (4) , and (8) governing 
the reactor kinetics is obtained by the same perturbation technique used 
to analyze the "T ^ model." This method investigates the effect of a 
small perturbation introduced by the feedback effect. The equilibrium 
state that ensues from the perturbation is defined by a series expansion 
about the initial unperturbed state of the reactor. The perturbation 
technique offers multiple solutions and unique eigenvalues for each of 
the equilibrium states (see Figure 2). The presence of more than one 
equilibrium state, known as the bifurcation phenomenon, requires an 
investigation of the stability of each of these solutions . The procedure 
of analyzing the stability entails imposing an initial deviation on the 
equilibrium flux and then examining the resulting temporal behavior of 
the deviation . 

A. EQUILIBRIUM STATES FOR THE SLAB REACTOR 

-3/2 

The first consideration to be made in the analysis of the "T 
model" is that of determining the possible equilibrium states of the 
reactor after a feedback perturbation has been imposed. The solution 
to the governing Eqs. (3), (4), and (8) is desired for the situation where 
the perturbed flux has achieved the equilibrium state (i.e., 2)) ^ “^)- 
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The governing equations reduce to the form 



= o 



O in £» <3) 

OH 



and 



dtp 



a 



; = lo 



( 8 ) 



The equation for the multiplication factor is easily shown to be 

'i(^) - §o , 

The function Vo constant are considered the 

initial flux distribution and multiplication factor before the feedback 
effect becomes significant. The result is substituted into Eq. (3), 
yielding thereby; 

2 r / \1 ( 71 ) 

V +L A + a “ ‘/lCx))J‘P(x) =0 ini> 

^(x) - O o^^ iD 

where 

A = L- 1. 

After the perturbation due to feedback the new equilibrium state 
is described by a Taylor series expansion about the initial flux distri- 
bution, ^ ^ , and the initial eigenvalues, o .If the amount 

of perturbation is measured by ^ (the perturbation parameter which 
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will subsequently be defined) , the resultant reactor state is 



^ (X)(z) — ^Cx) + ^ + 2 • (11) 

'X — Xo + €A X ‘ \ 

For small perturbations of the system the values of and Mo 

are described by truncating the series after the third term. The un- 
perturbed states for the reactor is defined as the solution to the problem 

V V^CX) + Xo ^o(x) = o in D (72) 

( X ) — O 6 2 > , 

It is readily seen that 

~ (2os V A. ^ "X 

/n = 1 , j, s, • • • , ^ 

So as to simplify the subsequent analysis the shape function / /h is 
defined for the slab reactor as 

= dos ^ . <”> 

Thus , 

f^Cx) - ^ . P4) 

The next term of the expansion is determined by differentiating 
Eq. (71) with respect to ^ and taking the limit as ^ tends to 
zero. This differentiation yields 
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= -(^(6) - f ‘?0‘,e) ‘4 ^cx^^6i) <fc\ 



t) It) J> 



^ C\b) - o 



Oy\ ci J) 



As ^ approaches zero, approaches t 0< j , XC6:) 

approaches A© / 

^^^C'X) ("Ao ~ z 

, in j) 

^ ( X) C> OK <iD, 

In order to obtain a solution to Eq . (76) the solvability condition (given 

in Appendix A) requires that the 



J r hf 

^ cx) cix =0 

T) 



(77) 



// 



where fCx) is the homogeneous solution to Eq. (76) given below 

z ■ // r 1 • (78) 

V + [■^^ “ 2 J =0 mi 

* // 

^ Cx\ - O OY\ . 

A solution to Eq. (78) is offered by a modification to the R. Courant 

and D. Hilbert perturbation method (provided in Appendix B) . Their 

approximated solution is given by 

* H (X r- /Ci\ ( 79 ) 

9^ CX) = U. (-X) -h Z 7^ (z) CX) -h • • • 
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Since the contribution of m is very small only the first two terms 

of the series will be retained. The substitution of Eq. (79) into Eq. (78) 

leads to the requirement that Un,(x) be the solution to 
2 . 

V U/v»C>C) + ^orA ^fA(X) — O 



UmOc) -= O 



OV\ 6 U 



the solution to Cx ) can be obtained as 



oo 









d 






A 






A 






J 7^ 



where 



ci /v\ 



= / 

Jy> 



-'/z 

(>{) U/A do< 



- 



The orthonormal solution for is LL/^Cx) 

Ao/w ~ ^ • The constant ^ /vxj 



A\ and 



determined as 



dMj = 



L A 



dx 



thus 



03 



1 



ol 



/v\ 



0- 






IS 



A^CxI =/_'ZT' 'J 

These results may be combined into the expansion given by Eq. (77) to 
yield the following result, 

CO 

a 



<pcl = 'V'j + 



AV 



(80) 
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For the convenience of the subsequent analysis this function will be 
referred to as the function given by 




A question may arise concerning the convergence of the series 
used in Eq . (81). If the Weierstrass M test for uniform convergence is 
employed, a series which is a continuous function of in the domain 

D is convergent under the condition 




where 






is a known convergent series of positive con- 



stants [6] . If the convergence of 




is doubted this test may be used. First, the maximum value of 
is determined as follows: 












d% 
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throughout the domain. 



Since the function 






/AJ. (j'-m) 



Ap 






If the convergence of the left hand side of the inequality is provided , 

then "Pj Cx) is convergent by virtue of the Weierstrass M test. The 

integral test provides a convenient method to establish this convergence, 

Hence, if the integral 
•00 



/ 

^ J 






exists, its convergence and that of ^ J is determined . 

AT - j J 

Thus, 






CO 



J 



_ 2_H\ 

Jl 



1 - 



where JT is any possible value of J . Since the integral exists, 
the series converges. 

A rs ^ 

The values of 9o ^ (X) from Eqs . (74) and (80), 

* 

respectively, are substituted into Eq. (77) and ^ is evaluated as 
shown. 
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( 



r 



Because of the orthogonality of the 
reduces to the following form. 



oo 



rZ ]dx=0. 

^ />v 






/VI function the above equation 







This equation may be satisfied by the condition 

^ = O. (82) 

The value of ^ is substituted into Eq . (76) and the homogeneous 

• H 

equation that results is similar to that defined by . Thus, 



ipiX) = 



S/» T^. 



(83) 



The value of the perturbation parameter, ^ , is chosen to represent 

the average magnitude of the linear term in the expansion for Eq. (11) 
weighted by the % function and measured in units of the maximum 
initial flux. Thus 

If Eq. (84) is differentiated with respect to ^ , the result is 

1 ^ ^ + ■ -Jdx . 
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By letting ^ approach zero and employing the solution for 
(i.e., ) the above equation reduces to 



^ _ r 

BA^ Jj. 



A» 



c/dc 



or, from equation (81) 



/q 



M 



B 



AX 




00 



-f 



2a aL. , 

j^L ‘j ^ 

j-W/- 

a* AX. 



] 






The quantity within the parenthesis is expanded and 



A 



zr» 



B 



/A 






A\ 









+ a/l|J' 



“ j . . nZ 

a Ir^ 



<L ‘j 



^•4=lj3,S, 



dx + 0 ). 



The second term on the right hand side is identically zero because of 
the orthogonality of the function. When the integration of the 

first and third terms is evaluated, i.e., x = 1 , the result 



IS 



Am 

P A\ 



h ^ 

— 1 + 0- 



oo 



l-Z 



dzAj. 



t2 



Tjv 



+ 



0(o?) 



^ c. 

Since the order of CX is generally very small, the contribution of 
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the higher order terms are neglected. Therefore, 



- Bm. ■ 

Equation (83) that defines ^ ("X) may be rewritten as 

<^(>0 = (85) 

The third terms in the expansions for and 

given in Eq. (11), are determined by differentiating Eq. (75) with respect 
to ^ and then taking the limit as ^ approaches zero. The dif- 
ferentiation yields 

/•• a •• -y^ 5a . \ 

As ^ approaches zero and Eq. (82) for is used the above 

equation simplifies to the form 
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fr 









I 9^X) 



= ( 86 ) 



in ji 



^ (y^) 



- O 



On ^ » 



The solvability condition for Eq. (86) is 

[ ^oCx) + ^ ^o<-y) ^cx) d'x — o 

^2) 



H 



(87) 



where Cp(x) » the homogeneous solution to Eq. (86), may be seen 
to have precisely the same shape function as the solution to Eq. (78). 



Thus , 



.. H 

^ Ck) 



n 



T 



( 88 ) 



• • • H 

If the values of . ^<rx) , and X) from Eqs . (74 ) , 

(85), and (88), respectively, are substituted into the solvability con- 
dition of Eq. (87), the following becomes evident. 







r /) '"V ^ 




The integrals are simplified using the definition of given in 

Eq. (81). Thus 
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f/y^ / 


^ r r 


\WryPy^ 


d X =/ '' 


L / u j 


K L 



oo 



'/>1 






If the contribution of any term having a coefficient of 



considered negligible, an approximate value of 
Hence , 



A 




Q. or larger is 
may be evaluated. 



^ C rA 






a. 









(89) 



The definite integral may be evaluated for the fundamental mode of 
reactor operation (i.e. , /V\- 1 ) and the value of for /Vl — i is 

\ ^ (. 279 ) . 

The equilibrium solution for and ^ (^) is obtained from 

the results of Eqs. (72), (74), (82), (85), and (90). The expansion for 
6) and ^ (^) , equation (11) is 

<P(%i) = (i'.. + € r„ i- 0(6^)) . 
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I 




I 



where 



are defined by Eqs. (73) and (81), respec- 



lively, and for the fundamental mode 

A (^) = 2 (.279) + ■■■ 

= iWz)^ + .01H9 + ..._( 92 ) 



The magnitude of is determined by the "power condition" 






0< 



where P is the reactor power level and ^ is a conversion factor. 

Equation (91) implies that there are two equilibrium states for each 
mode of operation corresponding to the negative and positive values of 
^ . For the fundamental mode this condition may be depicted as 

shown in Figure 2 . For reasonable values of ^ and ((X.) the con- 
tribution of Ql?/y\ is essentially that given by . The impact 

of the series of cosine terms within the Xa\ function may be seen to 



be negligible if the peak amplitude 
contributor is considered. 



£ f\ 



I 3 AX 



of the largest 



e A 



3 AW(X 



2_Q. 

J 



^ J-- 



d 



I 3 



1-3 



1 



X = O 



The value of ^ 13 ; ~ 



was calculated and for a very 

^ . 9^ 

conservative value of £ = . 1 and /S’, = -oox ■ 
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€/l 



13 M/qx 



= .1 



2 x.OOl -IIH4 1 



J 



1 8 ‘■ 



= 2.91 X iO 

« i . 



— 6 > 



For each of the equilibriurri states there are unique eigenvalues as shown 
in Figure 3. For a reactor with a negative coefficient of reactivity the 
fundamental eigenvalue is drawn as a parabola centered at (^/z) 
which opens to the right. Thus the equilibrium eigenvalues are always 
larger than those of the initial state. As has already been deduced for 
the "T ^ model" any decrease in the neutron flux (i.e. , negative ^ ) 

cannot be compatible with an increase in the material buckling. There- 
fore negative values of ^ for negative reactivity are unrealistic. 

For the case of positive reactivity the eigenvalues react in an opposite 
way as depicted by the dashed line. Values of ^ which are negative 
represent possible states of flux, but positive values of ^ are un- 



realistic . 



-1 



When the eigenvalues of the "T model" are compared with those 

-3/2 -1 

of the "T model" on Figure 3 it is noted that the "T model" is 

affected more by equal perturbations on the system. This result is con- 
sistent with that predicted by A. Reichel [3] as can be shown on Figure 1 
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There it is observed that the line for the "T model" is in very close 
proximity to Reichel's line of maximum Doppler effect. Hence, one may 

conclude that the Doppler effect has a more pronounced impact on the 

- 1 -3/2 

"T model" reactor than on the "T model" reactor. 



B . STABILITY ANALYSIS OF THE SLAB REACTOR 

The stability criterion for the equilibrium states of a slab reactor 
with negligible concentrations of delayed neutrons is developed as 
follows. The governing equations, Eqs . (3), (4), and (8) are combined 
into the equation 






= V 



“ O. ( 



(93) 



\v\ ]> 



•= O Oy^ 

^Ck,o) - ^cx;. 

The stability of an equilibrium state is investigated by prescribing 
an initial deviation, , in the equilibrium flux and observing the 

growth of this deviation. If the effect of the deviation increases the 
flux with time, the system is said to be unstable. If the effect is 
reversed the system is stable. 

The time dependent flux is described by a Taylor series expansion 
about the equilibrium state, , defined by Eq. (91): 
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( 94 ) 



such that 'l^^o . 

Equation (93) is differentiated with respect to and the result is 



u 


7 V, + [ 






]<f 


in ^ 




= O 




Ovi <£]> 








where 




7->J) . 










If the limit as tends to zero 


1 is taken, the equa 


tion reduces to 




Z>'t 




a 

z 


) - 


(X%U) 




(95) 

<n j) 




II 

0 




o>\ ^ D t 






The solution for 




is chosen to be 










k.'t) - u 


CK; t) 


[o( 




(96) 



Hence the stability of the system may be analyzed by inquiring whether 
o( (O is negative or positive. The value of oU^) is described 
by a Taylor series expansion about the parameter ^ and is given by 

• / 2L * ** 

^ oio ^ oL -i- ‘ ^ ( 97 ) 
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If the proposed solution for 6) (i.e. , Eq. (96) is substi- 

tuted into Eq. (95), then 

- q((^) ^ 

- CL 9i(x)J e) = O in i> 

C 6 ) — O o>» ^ i , 

Taking the limit as ^ approaches zero results in the expression: 

-.■2^ f . Ct- //5 \ (99) 

V U(7t^o) -f ( iXo - Z •= O <fj J> 

U.{X^C>) = O ov\ <^ ]> . 



A solution to this equation is offered by a modification to the R. Courant 
and D. Hilbert perturbation method (provided in Appendix B) . Their 
approximated solution is given by the series 



— Qr -7. , — (100) 

U. — UCOv) -f 1 N L'^) OJ~Cy.) ^ 

Since the contribution of is very small,only the first and 

second terms of the series are retained. The technique provides that 
U.CX) is the solution to 

) -h (A oon U-Cx; (An liA 



u t x; 



U ^ X) = o 



Or\ 
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Since the value of ^ has already been determined in Eq. (72), 

the only nontrivial solution for LA, ( X) is 






“• fYY\ / /v»0 

where is defined by Eq. (73). (Proof of this solution is pro- 
vided in Appendix C) . The eigenvalues corresponding to X) are 

A o <A OA^ . Therefore the first term for the series for 

0( U) is defined as 



o( 



OfY\ rm 



A A 






( 101 ) 



The linear term of the series of Eq . (100) is given by 



Co 



rvr 



-I 



frr\A — > 

— ^ LA. 



rf: /Vn 



\c>Dm Aoj 



An 



where 



c\ /w\ 









The solution for u may readily be seen to be 

^ A • 

tZ/ u. 



u = Ll 



Ta 



An I /An ' 2^ Ao/v»»'~Aoj 



7 






^ /yr\ 

d /vnj ~ ^ ^ ^ ' 



where 

This solution may be put into the form afforded by Eq . (81) 

LA. = ~LX. ^ /yy\ 



( 102 ) 
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If Eq. (101) is considered, some conclusions as to the stability 
of the reactor for very small € can be made. For the fundamental 



mode of the reactor ( /VI = ) the value of o/oi/m is 



oi on ~ ^ Oi 


^ O Avw 


• 


(l " /yvL ) . 


The possible values of olo\/yn. 




ol o\ ~ ^ 


/m = 1 




3/ 5 ; 7^ 



(102a) 



Thus for values of £ near zero the first term of the series for 

(^) indicates, that the flux is neutrally stable for the fundamental 

mode of the deviation \ ) and stable for the higher modes of the 

deviation (/WI?'l ). The question of neutral stability obliges an in- 
vestigation into the remaining terms of 

Some conclusive results may be obtained concerning the higher 
modes of reactor operation. The possible values of oU(=) for the case 
where An >1 are 










An 


> 






0 


/m 


= /n 


(102b) 




25 ^ 


An 


^ Al . 
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I 



I 



L 

) 



P 



Since o(ii) may be negative, the higher modes are unstable. 

The linear term in the series expansion of di(^) is evaluated 
by differentiating Eq. (98) with respect to €. and thus 

d. (K, t) + fe) - i Q. - a 

- uc-x,t) + [X(6) - - UCil 

. U ('X/^) = O . 

The value of has already been determined to be zero. If the 

limit as ^ approaches zero is taken, then 



+ [Ao ~ Z ^ Utx,0) (103) 

ucx,o) . 






The solvability condition requires that 

-Vz . 



/ ( a f\~ ' - \ w (i( 

( ^Cx) -t ixCX) Lx(x) dx = O 

« H 

where is the homogeneous solution to Eq. (103). If the value 

of ^^OAi/)n substituted into the homogeneous equation, the result 
is an equation similar to Eq . (99), thus 



(104) 



. H 

U 



H 



A. - 



When L-l />n and lA An, are used to evaluate the solvability 
condition the value of cK is developed. The details of this procedure 
are developed by using the results for ^^CX) . C K) , and 
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Li (X) obtained from Eqs . (74), (85), and (102), respectively 

. H 

and U (^) . Substituting these functions into Eq . (104) gives 



or ^ 

d'x T„Ty,„ d'x = o. < 1 “ 5 ) 

Consider the coefficient of C>( in view of the definition of 
given by Eq. (81) 





-f 

= i + O + OCa.’-) 

= 1 . 




d-x 

0(a’-) 



An approximation for the linear term of Eq. (105) may similarly be 



obtained to be; 




if terms of the second power of (a ) and larger are considered negligible. 
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These approximations are substituted into Eq. (105) and 

(X = ^ d-x . 



( 107 ) 



The stability of the system is determined from the evaluation of 
by substituting the results of Eqs . (101) and (107) into the expansion 
of Eq. (97). Hence, 







. (108) 
doc " 

» 



The additional term enables the further discussion of the stability 
of the fundamental mode of reactor operation. Thus, for /H = jL 




For all values of /YV\ 
the value of 0((€) 
defines the value of 



>X the first term becomes dominant and 
to be negative. However, the choice of 



(XCO , 



forces 

/TYl-i 



die) = 




(110) 



This integral may be evaluated and the result is, 

o((0 = -G ^ (.279) . (Ill) 



This result provides the following conclusions; 
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a) negative reactivity 



^ ^ for positive G 

> o for negative ^ 

These results indicate that the equilibrium solution described by 
an increase in the flux over the initial flux distribution will be stable. 
On the other hand, the other equilibrium solution, which was determined 
to be unacceptable, is unstable, 
b) positive reactivity 



0( (€.) 


o 

A 


for positive €. 


d (.O 


A 

0 


for negative C 



The equilibrium state described by a negative value of 6 
is observed to be stable. The other equilibrium state which was con- 
cluded to be unrealistic is unstable. 

In summary, the stability of a nonlinear reactor was studied by 
solving the boundary value problem and then considering the effect of 
an initial deviation on these solutions, called the equilibrium states. 
The possible equilibrium states are those described by Eq. (91) with a 
unique eigenvalue for each state given by Eq. (92). However, the 



d (€) 
d U) 
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stability criterion predicts that only one of these states are stable and 
that state will be physically observed in the reactor. For a negative 
temperature coefficient the stable state is one described by a positive 
^ . For the case of positive temperature coefficient the negative 

value of ^ represents the stable condition. Other equilibrium states 
are unrealistic. 



C. EQUILIBRIUM STATES FOR A FINITE CYLINDRICAL REACTOR 



The analysis of the equilibrium states of the finite cylindrical 
reactor can be done in a manner similar to that done for the slab reactor. 
The equilibrium states and the associated eigenvalues will correspond 
to those developed for the slab reactor. However differences will be 
noted because of the differences in the geometry. The flux distribution 
depends upon a radial coordinate, P , which varies from the center 
line to the outer radius, , and a longitudinal coordinate, ? , 

which varies from -1 to +1. 

The equilibrium solution is developed from the governing equations 



which are simplified by requiring 



c) 



2 )^ 



= O 



The governing equations, Eqs . (3) and (8), reduce to the form, 

2. r / ~~V'u \~i (112) 

iVv JO 



= o Ov» . 
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The perturbation solution is evaluated by a Taylor series expansion • 
about an initial reactor state, and . If ^ is the 

perturbation parameter, which will subsequently be defined, the resultant 
equilibrium state is 

< . (113) 



The initial states, , of the reqctor is defined as the solution 

to the system: 

= o in J) (114) 

= o on , 

The function is determined in a manner similar to that already 

provided for Eq. (43). The result is 



If the shape function 




is defined for the cylinder as 



% = dos('^5/z 2 ) 

then 

The eigenvalues are also determined by 

th 

where /v\ is the n root of the zeroth order Bessel function. 



(115) 

(116) 



( 117 ) 
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The linear terms of the expansion of Eq. (113) are determined by 
differentiating Eq. (112) with respect to ^ ■ and then taking the limit 
as ^ tends to zero. The result of these operations is 

+ (Ao - f ycos) = 



= o 



Ow (^i> , 



(119) 



In order that a solution is obtained for Eq. (118) the solvability con- 
dition (given in Appendix A) requires 

j 2jjr dr O 

where ^ ^^e homogeneous solution to Eq. (118). The 

.* // 

equation defining ^ 2 ) » Eq. (120), is solved by a modification 

to the R. Courant and D. Hilbert perturbation method (provided in 
Appendix B) . Thus, 

+(Ao - 2 ‘(iw) 



H 



//» 



^('• 7 ^) - o sy> 



The approximated solution is given by the series 

H a 



^cr,^) = U^^r.s) + 2 4- 0(Cl) .(121) 
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Only the first two terms of the series are retained since the contribution 

1 

of a is small. The solution for U^«7S) is obtained from the 
equation 

+ X, U„C'T,£) = o in T> 



= O 



OA , 



Since the solution for this equation was already determined for Eq. (114), 
the result is similar and 



Urn Cr,^i) 




The value of a /v\ is determined by the orthonormal condition im- 
posed on the 0./A O', function. Thus, 

2j r dr dz = i . 

The integral is evaluated and the result is that 

■ ( 122 ) 



The second term of the expansion is defined as 




^ Av. 
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where 



= j Uj(r,i)ljrdr di 



and 



Xo„ = (^-/ry * r%f . 



The value of ATm(^, 2) is 



Co 



/V,y,cr,^) 



-Z 






A = /,3,5/- 

^ /K. 






where 



d/Aj = aI^' 0/>n®j 2Tjj rdrdi. »23) 

. // 

The solution for ^ C^i^) may be determined by evaluating the 
series of Eq. (121). Thus, 



= B 



/H 









If the shape function J-AV is defined as 



T. . 



CO 



T„ = <S>^^ + ■ 024) 



/Vt 



v/here , A oA\ / A , 6nd are defined by 

Eqs. (122), (117), (123) , and (1 1 5) , the value of is 
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(125) 



= Sm Xw . 

The values of , given by Eq . (110), and Eq. (125) 

are used to evaluate the solvability condition , Eq. (119). Thus, 

BJ! '21, 2jdr dz - O 

or 

i I ^ dr dz = 0 . 

This equation is satisfied by the condition 

i = O . 026 ) 

« 

When this result is substituted into Eq. (118) the solution for 

is the homogeneous solution to Eq. (118) or the solution to Eq. (12 0). 

* 

Thus the value of (PCr,^) 

is 

^(r,^) = 6^ XL . (127) 

The value of the perturbation parameter, ^ , is chosen to represent 

the average magnitude of the linear term in the expansion of Eq. (113) 
weighted by the ^'/k function and measured in units of the maximum 
initial flux. Thus 

e = dK . 028 ) 
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r 

I 



If Eq. (128) is differentiated with respect to ^ and the limit is 



taken as ^ 

A 

Bav. 



approaches zero, the result is 



= f T 

/ D 

= /'r©A.'V'/i, + ®j 2jrdrdi 



/A 2tjt dr d^ 



i ^ ' 

/w 

If the quantity within the bracket is expanded and terms of ^ are 
considered small, then 

f\fA 



e 



A\ 



= 2rrr dr di 



+ a 






^ A^‘ 



+ 0 . 






Because of the orthogonality of the TA^ function the second term 
vanishes and 

Am 

B/y>‘ 

= 1. 



= y 0^ 2y r dr d2 



Thus, A /y\ “ 3/yx and Eq. (127) may be rewritten as 
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% 




I 



i 



( 129 ) 



The quadratic terms of the expansion for and 

given by Eq. (113) are evaluated by differentiating Eq. (112) twice with 
respect to £ and then let ^ tend to zero. The result is 

+(Ao ” t 

- ^ ^ Ao('Ci> A in 7> 

* • 

^ — O Ov\ <S1) . 



The solvability condition requires that 

Ol .n “Vz 



f / >• r. Ol ,ri~ ^ \ . ( 131 ) 

j ^ ^Cr,i)) (P(r,i:)l-^r d>r ~0 



H 



where Acr,i) is the homogeneous solution to Eq. (130). This 
was previously determined for Eq. (12 0) and from that result 

• t H ^ U ty\ 

(r,i) — C/n X/M. * 



(132) 



/ / 

The value of ^ may be determined by substituting the correct 
functions determined above . Thus, 
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X T/k " 

§ / Am ^’m. Tm 2 jr dr d?: 



The coefficient of ^ is evaluated as follows; 



Am / t 



2jr dr dz 



= Am/ tn f(9 + IZ >voJ- Xoj % t) 



A 



-X 2 j r dr d ^ 






U 






where a M is defined by Eq, (122) . The integral on the right hand 
side is evaluated as follows: 

^Am/'^m'^^^m 2jr dr di 

Am / <^M 2jr dr di 



rU Q- 



<L * * 

if terms multiplied by CX are considered small. The value of ^ 

A ^ ^IT'' dr di. 



(133) 
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For the fundamental mode (/n ) 



X = 



. S63 O. 









( 134 ) 



The various terms of the expansion of Eq. (113) may be combined 
to give the equilibrium states of the reactor. Thus, substituting Eqs . 
(116) and (129) into the expansion for gives 



6) = /)/M 



^ 0(i) 



(135) 



and substituting Eqs. (117), (126), and (133) into the expansion for 
A ) gives 




For the fundamental mode of the reactor (/n=l ) the integral is 
evaluated and 



>>('£)= (P'/rP) +(^z)^ + 6^ 



. S63 

ri- fA,' 



(136) 



It is apparent that for any mode of reactor operation there will be two 
equilibrium states defined by Eq. (135) and each of these states will have 
unique eigenvalues defined by Eq. (136). However, only one of these 
values will represent an acceptable reactor state. 
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D. STABILITY ANALYSIS OF THE FINITE CYLINDRICAL REACTOR 



The stability criterion for a finite cylindrical reactor with negligible 
fraction of delayed neutrons is determined by the methods previously 
considered. An initial deviation, ^ , due to temperature feedback 

is introduced into the equilibrium states of the reactor and the temporal 
effect is noted. 

The governing equation has the form given by Eq. (93) where the 
functional dependence of the variable is changed to accommodate the 
cylindrical geometry 



The deviation is described by a Taylor series expansion about the 
equilibrium states, Eq. (135), given by 




in D 



= O 



On ^ ^ 







such that 
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Eq. (13 7) is differentiated with respect to and then the limit is 

taken as the value of tends to zero. The result of these operations 



yields the equation 

“ I - 0-^0 



(139) 

in J> 



% O 



ip.. 



A form for the solution of is chosen to be 

^Ar,%i,'c) - u e/xp{_oU6)'c] 



(140) 



The stability of the system may be determined by investigating the sign 
of o(U) . If the sign of o/(^) were positive, then the exponential 
would increase with time. Hence, the contribution of in the 

series of Eq. (138) would be unbounded with time and the system would 
be unstable. On the other hand, negative value of would indicate 

that the contribution would die and stability would be assured. 

The value of (^) will be determined by a series expansion given by 



c((^) oL^ ^ Ck zi. U ^ 



(141) 



The following will be directed toward evaluating each of the terms of 
this expansion. 
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* 



If the proposed solution for % (r. , Eq. (140), is 
substituted intoEq. (139) the result is. 



+ [x^^) - o( (.a + f (1«) 

"'/i 1 

~ OL <Vo U. •= O ivi J) 

u. ( ^ ) “ O ov% <S J) 

where U. = U . The first term of the series of c^(i) 

is evaluated by taking the limit of Eq. (142) as ^ approaches zero. 

Thus , 

V’’u(Y:i; + [XonN - 

K u.c^ ^■) = O ivi i) 



U(/;^) = O ov% £P . 

The solution to this equation may be determined by the approximation 
given by the R. Courant and D. Hilbert perturbation technique. A 
solution to Eq. (143) is developed about the series 

4 j 4 - 0 i<k) . 

The solution for LA. C is determined about the solution to the 

equation 

a cr,i) 4 ( )\o,y. - O^o) - O m3) 

d cr,^) =0 on sr> . 
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It may be noted that the solution to this equation is similar to that ob- 
served in Eq. (114) . However, by the results of Appendix C one may 



conclude that the solution for 



Lx (f] is 



where the function 
Thus , 




is defined in 



given by 

Eq. (115) and X O/Wl ^OM O(o- 







c rr\ 



X 



o/vr\ 



(144) 



and the linear term of the c<(i} series of Eq. (141) is defined. (A 
discussion of o( e> will subsequently be made .) Since the Lx(«72) 

function must be orthonormal the constant u />n rnay be evaluated 
by the normality condition 



L 






/ir^ ' /m 



Ijr dr di: = i . 



The value of u /yy\ is 






AW 



0 



z 

AA. 



TJ ri 



(122) 



which has previously been determined by Eq . (122). Thus LX. "i) 



IS 






The second term of the expansion UXC'T, S) is defined as 
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OO 

Z d AV\ -1 — 



or 



^ />w. 



CO 



«->■:=> “I -ftf’-., 0)'^) 



i = '.5/S‘/' 
4 /r^ 



where cImj is defined by Eq. (123). The solution for the 
function from Eq. (143) is determined as 



U(^^) = 


UCy;?:) 


C — 

2 ^ ' • 




L- 


oa 

£ y d 

2 - ^ OAA ~ <^0 4 


= T-X /vn 


L 



] 



(145) 



where ^/w\ is defined by Eq. (124). 

The value of oL o/«/r». as evaluated from Eq. (144) is rewritten 

below 

0(or*\^ ^ O AA ^ O/WV. * 

If the value of ^ 

For the fundamental mode ( M - 1 ) this equation reduces to the form 



(146) 



o as given by Eq. (117) are used, the result is 

i / 7 -1 \ I .»T-. \7- 



(146a) 
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o(o,^ = at (^& - 



For values of £ near zero the values of 



/7TL are 



oC oi 



/V*y 



= 0 
< -Z*j 



/yn - /yK. 
/m > 



Thus for all modes of the deviation except the fundamental the sign of 
0( OAl is negative and any initial deviation of the system will 
decay with time. However, for the value /W\-A\ the first term evalua 
tion of the stability indicates that there is a neutral stability for the 
fundamental mode. Thus it is necessary to investigate the higher order 
terms of the expansion of . For the higher modes of reactor 

operation o(U) may have the following values, 

d CO ; 



25' 


/yri ^ /K 


0 ' 


/m ■= AX, 




Ayx. > AX- 



Thus the value of dU) may be positive when M for //i. > 1 

Since the value of />v\ Is so large when compared to the linear 

term (order of £ ) of d U) , one may conclude that the first 

term of the series is dominant and the higher modes are unstable . 
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The linear term in the expansion of d(a , Eq. (141), may 
be evaluated by differentiating Eq. (142) with respect to ^ and taking 
the limit as ^ tends to zero. Thus, 







I / 1 . 

M O “ 2 ^ To ^OAl/W» J U. (V, 2) 

“Vi * , \ 

d + k) 



in 3) 



u cr,^) “ o c»w ip , 



The solvability condition requires that 

^ / -J. * \ ^ 

/ (a + ^ ^ uCo^) 2jrdrd^ 

. // 

where UL is the homogeneous solution to Eq . (147). This 

solution was previously developed for Eq . (143) and the solution is 



= o 



UiCr,^) = ^/m • 

. // 

With the solution for ^'72J and the foiutions for , 

# 

, and given by Eqs. (145), (116), and (129) 

the value of (A may be evaluated from the solvability condition. 
Thus , 
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(X U /m ah X m 2ijr dr dz 

- ^ f\!^ ^ 'Xa^ Tat. Ijrdrdi 

The coefficient of 6 ( was evaluated for Eq. (128). The result in- 

z 

dicated that if the terms of order a are considered very small, 
the integral was equal to one. The value of o( is 

(i" = ■? 2jr dr di . d«) 

Since the fundamental mode of the reactor is of greatest interest, the 
integral of Eq . (148) is evaluated for /WT =r yvx. — j, . Thus, 

oil = "^ T/ 2jr dr dz 

r© 

The results of Eqs . (146) and (148) may be combined into the expansion 
of ciLi) ,Eq. (141), 
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( 149 ) 



oL 




dr dl 



and for the fundamental mode 

21 f 2L 

(O = (r;) \2.H0'^S - +(^z)(l-/>n) 

^ 3.'73 a ^ / i\ 

" ^ r; + 0 <:£ ) . (150) 

From the previous analysis values of /V>\ > 1 have been found to 
favor stability. For the value of /W = j. the following becomes 
evident, 

0( ^ ^ + OU) . (‘ 51 ) 

For the case of negative reactivity the sign of d ((:) will be deter- 
mined by the value of ^ .A positive value of ^ will permit a 
negative sign of c2(^) • The deviation associated with o/ (^) 

will decrease with time and stability results. For positive reactivity 
the sign of ( a ) is changed and the negative value of permits 

stability. The other states of reactor equilibrium are unrealistic as 
determined by previous arguments . 



95 



ippm 

ir 






APPENDIX A 



THE SOLVABILITY CONDITION FOR NONHOMOGENEOUS 
BOUNDARY VALUE PROBLEMS 

Given the equation 

^ - fcx) iy> 1> 



^(rc) = O 



OV\ 






prove that this equation has a solution if and only if the function 
satisfies the condition: 

’f('x) dV = O 






where 



1 r\ ^ ^ 

V t(^) - O in 7) (a-2) 



H 



^(x) ^ O 



Ovn 






Multiplying Eq. (A-1) by ^(X) and integrating throughout the 
domain the following results: 

V ^(x) + y)^(x) dV 

= f ^(x) f(x) dV , 



(A-3) 
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From Eq. (A-2) V = ->v and this is substituted 

into Eq. (A-3) . Thus, 

- <Pix) dV 

= f ¥’(x) -P(x) dV . 




By Green's Theorem the left hand side of Eq. (A-4) may be transformed 
into a surface integral; 



n -cfiu) v^%)) dS 

“ / 'f(x) ^(x) dV . 

Thus, if <p(x) is to be the solution to Eq. (A-1), the homogeneous 
boundary condition on S)> must be satisfied. Therefore, 



fcx) ^cx) dV - O 



(A-5) 



Q.E.D. 
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APPENDIX B 



A PERTURBATION METHOD FOR 
LINEAR EIGENVALUE PROBLEMS [7] 

The solutions to equations of the form 

Z _ - 

V U^v - ^ = O ii^ T) 

Ua\ = O oa iP , 

where € Y"(X) represents a small perturbation^ (|^ V'C'x)| ^ X a\ ) / 
are provided by a series expansion about the solution of the equation 

X/K U /K = O (B- 2 ) 

LA /v\ “* (S d • 

The eigenvalues are held constant while the eigenfunctions are perturbed 
in accordance with the Taylor series given by 

Um = U/v, + ^ + • • • _ (B-3) 

If the series is substituted into Eq. (B-1), the result is 

^ { U. /A ^ ^ 

Xa\ (U/y\ + CaT’a^ + LX /a 

- o . 



98 





r 



By separating terms of the same power in 



^ , the following equations 



are obtained: 





LA /y\ 


+ 


^/w 


La /v\ 


o 




(B-4) 


v" 


rsTZ 


4- 


^ /n 


- 


Yi'x) 


LA m 


(B-5) 




ur^ 


4- 


/V\ 


rc'K) 


LA/y\ , 


(B-6) 



Multiply Eg. (B-5) by LA and integrate over the volume, dV : 

y W rr^ d\/ r(x)u.^ixji dv 



j V y it AT/v» y U.ji^ dV 



(B-7) 






i 



rw)u,v.u^ dv. 



By Green's Theorem the following volume integral may be transformed 
into a surface integral 



/ (vVm “ y^ujjj dv = 

Jb ^ 



I n (UiVrr^ - /^r« Vujj dS . 



However, the boundary conditions specify that Li jp and A^v 
are everywhere zero on the surface and the surface integral vanishes. 
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Equation (B-7) may be rewritten as 




+ d V 



where 

cIms = f dv. 

Equation (B-2) reveals 

V Ujj = - \ Uji, 



d/v.j? 



(B-8) 



which is substituted into Eq. (B-8) to give 

()\/y\ — ~~ cImR 



where 



CX fr\ X 



j /^/y\ IX. jt 



dV 



(B-9) 



The value of the integral 
normalizing condition 
clear that 




Q /«/K 

Z 



is evaluated if the 
u! dv = 1 is imposed . It is 




c An 



An 



1 

Z 






-f 



dv - 1 



(B-10) 



By equating powers of ^ on the right and left sides of Eq. (B-10) 

J d'x i 

J /w d’X " o 



100 



or 



6 = o . 



The function is determined by a series of orthogonal 

functions given by 



oo 






(B-11) 



Multiplying both sides by and integrating over the domain 



/; 



Uj( dv 






. 



Assuming that the expansion exists and that the integral of the infinite 
series is equal to the sum of the integrals, 



CO 



CImA “ X /jN 



C[ AA A\ * 



Since 



and 



uy = o 

/ Uj dv = 1 









the constants evaluated as t> ~ OL/v\>\ . 



;3 
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However, from Eq. (B-9) it is obvious that K 



d M ^ 



Thus the series of orthogonal function given by Eq . (B-11) may be 
evaluated as 



oo 






d 



:)■ 






^ A\ *■ ^ J 



UxCx) . 



The series of Eq. (B-3) may be evaluated as 



05 



Z _d^J 



■ UaC-x) -t- 



j=17.y. 



where 






I 
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APPENDIX C 



Given the equation 

+ >>/v, 'Vai = o in ■& (C-1) 

Of ) = O on i x> 



with eigenfunctions eigenvalues ^/v\ • Prove 

that if 

-=0 ini) (C-2) 



^(X) = O 



OVk 4"I) 5 



where A >i /v\ is a new physical parameter, then the nontrivial 
solutions to Eq. (C-2) are 



where /y>7 and /Vl are independent subscripts. 

Multiply Eq. (C-2) by CX) and integrate throughout 

the domain. 

^ ^ (X) vVfx) dv + dV = o 



103 



(C-3) 



By adding and subtracting ^(x) i'^i.x) 

Cx) V^(fcx) - ^cx) V^^(x)) dV 

^ ^(x) dV +xj 't,(x) <P(X) dV = O . 

From Green's Theorem the first integral may be transformed into a surface 



integral by 




^(x) V ‘P(X) - ^(-k) dV 

= f n V ^(x) - cfcx) V fx;) d S . 

h 

However, from the homogeneous boundary conditions of 
on the surface integral vanishes. Thus Eg. (C-3) is 

I<P(X) (X) d V + > f Cx) <P>M dV = O r'’ 



From Eg. (C-1) 
therefore 



(X) = - 'V'mCXJ ; 



or 



- ii,(x)dv -I- (X) ^cx) dv = o 

(x - ^/«) / 



Cx) dv - O . 



If (X- \ ) o , then the integral must vanish throughout 

the domain and ‘P(x) is orthogonal to CX) . Hence, 
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^(x) — ^AY) 



and th6 eigenvalues of Eg. (C'*2) (i.e. , ^ 



>> = A 



/yyy 



Q . E . D . 



) are 
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